// ---------------------------------------------------------
// Class: MPC_Base
// MPC base class, modified based on Gangshi's code
// Assume the following feedback information is available
// - ht
// - rough2 (roughness)
// - g_dens (global density)
// - alpha and beta
//
// MPC_Base provide the necessary I/O functionalities. It does not set up particular 
//   optimization problem. In current version, MPC_Base comes with a optimizor 
//   (fminsearch).
// ---------------------------------------------------------
#ifndef MPC_BASE_H
#define MPC_BASE_H

#include <fstream>
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <iomanip>
#include <string>
using namespace std;
#include "../../GasPhase/RateLaw/RateLaw_Kleijn1989.h"
#include "../../KMC/plant_1DTRI_RDM.h"
#include "../model.h"

// ---------------------------------------------------------
const double TOLERANCE	= 1e-15;
const double VERSION	= 1.2;
#define PI				(3.141592653589793)
#define MODE			(100)
#define	DEP_HIGH		(0.45)			// add safety space of 0.00001
#define	DEP_LOW			(0.10+0.000001)	// add safety space of 0.00001
#define T_HIGH          (1000)
#define T_LOW           (300)
// ---------------------------------------------------------
struct MPC_diagnostic
{
	double  cost;
	double* Hopt;     // Predicted H, layer
	double* Dopt;     // Predicted SOR
	double* R2opt;    // Predicted R2
	double* varR2opt; // Predicted varR2
	double* M2opt;
	int     exitflag;
	double* ub;		  // Upper bound
	double* lb;       // Lower bound
	double  M2_0;
};
// ---------------------------------------------------------
class MPC_settings
{
	double Dens_set;
	double Fact_D;
	double Rough2_set;
	double Fact_R2;
	double Ht_set;
	double Fact_H;
	double VarR2_set;
	double Fact_VarR2;
	double M2_set;
	double M2_fact;
	double Rel_err;
	int    Horizon;
	double *Weighting;
	int    ControllerID;
public:
	MPC_settings(int iHorizon);
	~MPC_settings();
	int set(char* VarName,double VarValue);
	int set(char* VarName,int VarValue);
	int set(char* VarName,char* VarValue);
//	double get(char* VarName);
//	int    get(char* VarName);
//	char*  get(char* VarName);
	void PassToMatlab();
};


// ---------------------------------------------------------
class MPC_Base
{
protected:
	// SETTINGS
	double	t_int_ctrl;	 // sampling interval, s
	double	t_begin;	 // s
	double	t_end;		 // s
	int		Step_num;	 // Prediction horizon
	double	Dens_set;
	double	Fact_D;
	double	Rough2_set;	 // layer^2
	double	Fact_R2;
	double	Ht_set;		 // layer
	double	Fact_H;
	double  M2_set;
	double  M2_fact;
	double  Fact_varR2;
	int		Rel_err;
	double	T_ini;		// K
	int		Detail;
	int		ControllerTypeID;
	char	*FileName;
	char	FileName2[256];
	double  *uplim;		 // layer/s
	double  *lolim;		 // layer/s
	double	RT_T;		 // layer/s^2
	char    *Optimizor;

	// FEEDBACK INPUTS
	double	*alphas;
	double	*betas;
	double	ht;			// layer
	//double	rough;	// layer (NOTE: this is NOT r^2)
	double  R2;
	double	g_dens;
	double	P;			// Pa
		
	int		Repeated;
	int		t_idx;
	int		t_node;
	double  *t_arr;		 // s
	double	*W_sum;	 	 // layer/s
	double	*W2_sum;	 // (layer/s)^2
	double  *T_sum;
	double  *T2_sum;
	double  *v_dep_rate; // layer/s
	double  *v_T;        // K

	struct MPC_diagnostic MPC_diag;
	double cost;

	// MODEL:
	model  *Model;
	ofstream fd_time;
	ofstream fd_T;
	ofstream fd_W;
	ofstream fd_cost;
	ofstream fd_Hopt;
	ofstream fd_Dopt;
	ofstream fd_R2opt;
	ofstream fd_M2opt;
	
	// Some temperary variable for debug
	double   db_g_dens_SOR;
	void    log_detail();

	// DEFAULT OPTIMZATION SOLVER
	int		amoeba_constr(	double **p, 
							double y[], 
							int ndim, 
							double ftol,
							double para[], 
							int *nfunk,
							double uplim[], 
							double lolim[]);
	double	amotry(	double **p, 
					double y[], 
					double psum[], 
					int ndim,
					double para[], 
					int ihi, 
					double fac);
	double	fminsearch(	double p0[], 
						int ndim, 
						double ftol,
						double para[], 
						int *nfunk, 
						double uplim[], 
						double lolim[]);

	virtual double	cost_funct(double DEP_RATE[],double para[]);
	virtual void	upperlimit();
	virtual void	lowerlimit();

public:
	// Output
	double	dep_rate;	// layer/s
	double  T;			// K
	double	time;		// second

	// 2009.08.01 Gangshi: Add iDetail to regulate the output, 0: no detail files; 1: detail files
	// 2009.08.05 Xinyu: Add iControllerTypeID to indicate what controller is used
	MPC_Base( double	iDens_set,
			  double	iFact_D, 
			  double	iRough2_set,
			  double	iFact_R2, 
			  double	iHt_set,	
			  double	iFact_H, 
			  double	it_end,
			  double	iDepRate_ini = 0.1,
			  double	iRT_Dep = 0.1,
			  double	it_int_ctrl=1.0, 
			  int		iStep_num=10,
			  int		iRel_err=1,
			  char*		iFileName = "MPC_Base.dat",
			  int		iDetail = 0,
			  int		iControllerTypeID = 1);
	virtual ~MPC_Base();
	//double	update(	double ialphas[], 
	//				double ibetas[], 
	//				const double iG_Dens,
	//				const double iRough, 
	//				const double iHt,
	//				const double P,
	//				const double T);
	// Gangshi: use structure lattice_state to transfer state
	virtual double	update(	struct lattice_state iLState,
					        const double P,
					        const double T) = 0;
	virtual double  update2(const double time,const struct lattice_state iLState) = 0;
	virtual void	reset();
	int		        notfinish();

	inline double   getTime(){return time;}
	inline double   getDepRate(){return dep_rate;}
	inline double   getT(){return T;}
};
#endif